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ABSTRACT 

We numerically study the mutual interaction between dark matter (DM) and Popu- 
lation III (Pop III) stellar systems in order to explore the possibility of Pop III dark 
stars within this physical scenario. We perform a eosmological simulation, initialized 
at z 100, which follows the evolution of gas and DM. We analyze the formation 
of the first minihalo at z ~ 20 and the subsequent collapse of the gas to densities of 
10^^ cm~^. We then use this simulation to initialize a set of smaller-scale 'cut-out' 
simulations in which we further refine the DM to have spatial resolution similar to 
that of the gas. We test multiple DM density profiles, and we employ the sink particle 
method to represent the accreting star- forming region. We find that, for a range of 
DM configurations, the motion of the Pop III star-disk system serves to separate the 
positions of the protostars with respect to the DM density peak, such that there is 
insufficient DM to influence the formation and evolution of the protostars for more 
than ^ 5000 years. In addition, the star-disk system causes gravitational scattering of 
the central DM to lower densities, further decreasing the influence of DM over time. 
Any DM-powered phase of Pop III stars will thus be very short-lived for the typical 
multiple system, and DM will not serve to significantly prolong the life of Pop III 
stars. 
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1 INTRODUCTION 

The formation of the first stars was a crucial point in the 
early Universe. Believe d to have formed at z > 20 within 
lO^Mp) minihaloes (e.g. Haiman et al. l ll996l : lTegmark et all 
ll997l : iYoshida et all [20031 ') . these Population III (Pop III) 
stars began the process of transforming the Universe from 
a neutral and metal-free state to t he ionized and metal- 
enriched state wo observe t oday (e.g.|B arkana fc Loebl | 2Q0i ; 
Bromm fc La,rson 2 04: Cia rdi fc Ferra ra 2005: lGloveill2005 ; 
Bromm et a'iTT2009l : iLoebl |2010| ~ However, the extent to 



which the first stars drove the evolution of the early Universe 
is highly dependent on their mass, luminosity, and effective 
temperature. While the eventual mass of the first stars has 
usually been assumed to be driven by the cooling and chem- 
istry of the baryons within early minihaloes, recent work 
has posed that dark matter (DM) may also affect the mass 
and evolution of Pop III stars beyond simply providing the 
initial gravitational potential w ell in which the star-f o rming 
gas first collapses. For instance. ISpolvar et all (|2008l ): lloccd 
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(|2008l ) find that, if DM is composed of weakly interacting 
massive particles (WIMPs), energy from DM annihilation 
may also play an important role. 

DM annihilation may first become significant during 
the initial collapse of gas within a minihalo. This is due 
to the growing density of DM in the center of the mini- 
halo as it responds to the growing potential well of the gas, 
gene rally termed 'gravit ational contraction .' Several studies 
Ce.g.lSpolvar et al.ll2003 : iFreese et aT]|2008l : iNataraian et all 
I2OO9I ) find that gravitational contraction leads to sufficient 
DM annihilation to halt the collapse of the primordial cloud 
before a hydrostatic object has formed, leading instead to 
the formation of what has been termed a 'dark star,' a gi- 
ant (~ 1 AU) star powered by DM annihilation instead of 
nuclear burning. 

Because of the extended nature of these objects, their 
effective temperatures are too low to emit ionizing radia- 
tion. Depending on how long gravitational accretion of DM 
continues, this may allow for a much longer gas accretion 
period before the dark star phase ends, after which the star 
begins contraction to the main sequence, and radiative feed- 
back shuts off mass inflow onto the star (e.g. ISpolvar et al.l 
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I2OO9I ). Ilocco et al] (|2008l '). find that the dark star phase is 
short-lived, though they did not follow the gradual accre- 
tion of gas onto the s tars over time. More detailed work 
iRipamonti et al]|2010l ) even finds that such 'gravitational 
accretion' of DM during the initial cloud contraction will 
not halt the gas collapse, so a dark star of this type does 
not form. Thus, the initial mass of the first stars will not be 
significantly affe cted. However, oth er work reaches different 
conclusions (e.g. iFreese et al1l2008l ). 

At later stages. Pop III protostars may continue gath- 
ering DM through continued gravitational accretion as the 
protostar's mass and potential well grows. They may also 
gain DM through 'scattering accretion,' a process in which 
WIMPs scatter off the gas of the star and become grav- 
itationally bound to it. If the scattering cross section be- 
tween WIMPs and baryons were large enough, and the re- 
sulting capture rate of DM by Pop III stars sufficiently 
high, this would prolong the lifetimes of Pop III stars. 
This is because hydrogen will burn at a reduced rate while 
DM annihilation helps to su pport the star. Recent work by 
ISivertsson fc Gondolol (|201ll '). however, find that this phase 
of scattering accretion will be very short-lived, < 10^ yr, 
much smaller than the lifetime of the star. 

However, the extent to which DM will alter the nature 
of Pop HI stars depends highly on the DM density profile 
and the location of the DM density cusp with respect to the 
Pop HI star Previous work has found that a sin gle Pop HI 
star will form at the cen t er of a minihalo (e. g. lAbel et akl 
l2002l : iBromm et aLll2002l : lYoshida et "al]l2006l ). This is the 
picture that has been used by the previous studies of dark 
stars, and DM effects are likely to be strongest in this sce- 
nario. However, recent work has found that Pop III stars do 
not necessarily form in isolation. Instead of a single stellar 
peak at the center of a minihalo, a stellar multiple system 
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et al . 2011b). In particular, 
(|2010l) . which was started us- 
ing cosmological initial conditions, found the formation of a 
Pop III multiple system within a rotating disk structure of 
order 1000 AU in radius. Recent calculations have further- 
more found disk fragmentation and multiplicity to be robust 
against both photo- dissociating and io nizing feedback from 
the protostars (e.g. ISmith et al]l201ll : Stacy et al. 2011b). 
Even when allowing for stellar mergers on unresolved scales, 
Stacy et al. (2011b) find that a massive binary remains at 
the end of their 5000 yr calculation. This modified picture 
of Pop HI star formation is als o robust against statistical 
variation of the host minihaloes (jCreif et al.ll20lil ). 

Only one-dimensional simulations have been used thus 
far to study the possible effects of DM annihilation on the 
first stars. Three-dimensional effects, in particular stellar 
multiplicity and motion of the disk gas, may have impor- 
tant consequences for how DM and Pop III stars in teract 
in th e center of minihaloes (see discussion in, e.g., Iloccd 
One may initially guess that a DM profile could 
not remain peaked under the gravitational scattering ef- 
fects of a star-disk system, and would thus be less prone 
to scattering accretion onto the stars. As discussed in, e.g., 
ISivertsson fc Gondolol (|201ll ). however, the motion of a star 
through a DM halo may enhance its overall accretion rate. 



particularly for a lower-mass star that would need less DM 
to be powered as a dark star. 

To address these open issues, we perform a three- 
dimensional simulation which follows the evolution of a 
range of peaked DM density profiles in a minihalo that 
hosts a Pop III multiple system, yielding more accurate con- 
straints on the influence of DM annihilation on the nature 
of the flrst stars, as well as the influence of a Pop III star- 
forming disk on the inner DM proflle. Representing each Pop 
III star that forms with a sink particle, we follow the gas and 
DM evolution for 2 x 10** yr (~ 500 free-fall times for our sink 
density n = 10^^ cm~^), and we record the density of DM 
within the sinks as the evolution of the gas and DM profile 
is followed. We initialize our simulation such that the maxi- 
mum DM density is in the region of the most massive sink, 
giving an upper limit to the WIMP annihilation heating rate 
in this region as the star grows, as well as the potential for 
DM capture to prolong the life of the Pop III star. 

In Section 2 we describe the initialization of our simula- 
tions. In Section 3 we present our results, including descrip- 
tions of the DM evolution, the Pop HI growth rates, and 
the expected effects of DM on primordial gas and stellar 
evolution. We conclude in Section 4. 



2 NUMERICAL SET-UP 

The numerical simulations were run with GADGET2, a 
widely-tested three-dimensional N -body and SPH code 
(|Springel et al.ll200ll : [Springelll2005h . We used adaptive grav- 
itational softening for both gas and DM particles, but with 
an imposed minimum softening length (see Section 2.2). The 
gravitational softening length was thus assigned to each par- 
ticle based upon its current density, such that the softening 
length was adapted both throughout the simulation box and 
over time. 

A scale over which gravitational forces are 'softened' 
is necessary when representing a smooth mass distribu- 
tion with a finite number of simulation particles. The ideal 
softening length e to use in N-body and SPH sim ulations 
has been studied by a number of authors (e.s. [Merrittl 



1 19961 : iBate fc Burkertll 19971 : lAthanassoula et al.il200of r If e is 

too small, the gravitational force between discrete N-body 
particles can become arbitrarily large, and computational 
timesteps prohibitively small. Unphysically large fiuctua- 
tions in force among the finite number of N-body particles 
will result. If e is too large, real features on scales smaller 
than e will not be resolved (|Merrittlll99"^ ). For SPH parti- 
cles, collapse of Jeans-unstable clumps will be inhibited if 
e is greater than the smoothing length. On the other hand, 
setting e to be less than the smoothing length may result in 
artificial fragmentation. In a simulation with a large range of 
densities, the ideal e will not be the same in both dense and 
diffuse regions. Varying e with the smoothing length of each 
particle provides a way to avoid the problems of e values 
which are too large or small (e.g. iPrice fc Monaghanll2007l : 
llannuzzi fc Doladl201iri . This does lead to some extra com- 
putational cost in our simulations, however, as full adaptive 
softening requires smoothing kernels to be calculated not 
only for the SPH particles but also the DM. 

Our simulations evolve the central regions of the first 
minihalo that formed in a previous, larger-scale cosmological 
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simulation. We furthermore replace the DM from the mini- 
halo with more highly refined particles, assuming a range of 
density profiles. The gas from this region is al ready resolved 



down t o scales of 50 AU. As found in, e.g., Spolvar et al.l 



l|2008h : lNataraian et al] (|2009l '): lRipamonti et all (|2010l ') 



resolution length is the approximate scale at which the DM 
annihilation heating rate may first surpass the gas cooling 
rate and prevent gas collapse to stellar densities, given a 
sufficiently peaked DM profile that is aligned with a sin- 
gle stationary star. Our simulations thus have sufficient gas 
and DM resolution to determine whether this may still be 
the case within a Pop 111 star-disk system. In the following 
sections we describe in further detail our numerical initial- 
ization. 



2.1 Gas Initialization 

The gas distribution and its thermal and chemical state were 
taken from a sn apshot from a prev ious simulation of Pop III 
star formation (jStacv et al.ll201Cl ). This earlier simulation 
was initialized at 2: = 99 in a 100 kpc (comoving) peri- 
odic box. This was done in accordance with a ACDM cos- 
mology with = 0.7, Qm = 0.3, Qb = 0.04, and h = 0.7. 
To accelerate structure formation, we used an artificially en- 
hanced normalization of the power spectrum, ag = 1.4. 

This previous simulation employed a standard hierar- 
chical zoom-in procedure to i ncrease the resolu tion of the 
central 20 comoving kpc (see IStacv et al.l I2OIOI for further 
details). The most highly resolved DM particles had masses 
of muM = 0.096 Mq. Thus, scales of < 10* AU which en- 
closed < 1 M0 of DM (~ ten DM particles) were not well- 
resolved (see Fig. [1]). The most highly resolved gas particles 
had a mass mspH ~ 0.015 M0, giving a mass resolution 
of Mroa — 1.5A'ncighn^SPH < 1 Mq , where A''ncigh — 32 is 
the typi cal number of particle s in the SPH smoothing ker- 
nel (e.g. lBate fc Burkert|[l997l '). The gas evolution was then 
followed until it reached a maximum density of 10^^ cm"'', 
at which point we implemented the 'sink particle' method. 
This is a computational technique in which a small, high- 
density, gravitationally collapsing region is replaced with a 
single sink particle that will grow in mass as gas particles 
continue to fall onto this region. This circumvents the need 
to continue following the evolution of the star-forming region 
to yet higher densities, which is very computationally expen- 
sive, and allows for a numerical representation of a growing 
Pop III star. Our particular impleme ntation of the sink par- 
ticle method is identical to that in IStacv et all (|201ll ). In 
short, we check each gas particle to determine if it satisfies 
the criteria d < race and jsph < Jcont, where d is the dis- 
tance of the particle relative to the sink, race — 50 AU is the 
accretion radius, jsph = Wrotd is the angular momentum of 
the gas particle, Vrot is the gas particle's rotational velocity, 
Jeent ~ \/GMsink''aee IS the augular momcutum required for 
centrifugal support, and Mgink is the mass of the sink parti- 
cle. If these criteria are satisfied, the gas particle is removed 
from the simulation, and the mass of the accreted particle is 
added to that of the sink. Our sink accretion algorithm also 
allows for the merging of two sink particles if the smaller 
sink satisfies the criteria listed above. 

The gas of the central star forming region is taken from 
this previous simulation at the point immediately after the 
first sink has formed. The gas density profile at this point can 



be seen in Figure[T] Depending on the particular DM density 
profile, we cut out a central cube that is either 1 or 20 pc on a 
side, equivalent to a total of 3.6x10* gas particles (540 Mq) 
and 2.5x10'' gas particles (4000 Mq), respectively. Note our 
cut-out regions contain only the most highly-resolved gas 
particles from the original cosmological simulation, so every 
gas particle has a mass of 0.015 Mq. We then resimulate only 
the cut-out region, excluding the outer less-resolved regions 
that were in the full-scale cosmological simulation. The new 
cut-out simulations have the same mass resolution for the 
gas (< 1 Mq) and use the same sink particle technique, 
but now include the addition of a peaked DM profile that 
is resolved on sink scales and centered on the most massive 
sink. 

Note that the gas in the outer parts of the cosmological 
simulation can be safely ignored for our purposes due to the 
long dynamical times. The gas at the edge of the cut-out has 
typical densities of 10^-10* cm~^, corresponding to free-fall 
times of ~ 5 X 10^ to ~ 5 x 10*^ yr. This is over ten times 
longer than the length of time over which we evolve the gas 
in the cut-outs (20,000 yr). Furthermore, a rarefaction wave 
at the cut-out edge could possibly develop due to the vacuum 
boundary conditions. However, such a wave will only travel 
a distance of Cs t ~ 8500 AU, where Cs is the gas soundspeed 
(~ 2 km s~^), and the time t is 20,000 yr. This is less than 
one-tenth of the half-length of our smallest cut-out. 

By adding the highly resolved DM only after the first 
sink has formed, we are assuming that DM annihilation did 
not have a significant effect on the gas evolution at low densi- 
ties and was unable to halt the gas collapse before it reached 
densities of 10^^ cm~^. This is consistent with the results of 
iRipamonti et al] (|2010l ). who find that up to n ~ lO'^* cm""' 
the evolution of temperature with density will not vary sig- 
nificantly for a range of DM parameters. They furthermore 
find that heating through DM annihilation surpasses cool- 
ing at critical densities of ricrit ~ 1 0^ — 10^'^ cm~ ^ , sim i- 
lar to the critica l dens ities found bvlSpolvar et al] (|2008|) : 
iNataraian et al] (|2009l ). However, IRipamonti et al] (|2010l ) 
find that this will not halt the continued collapse of the gas, 
with the exception of a very brief (~ 3yr) stall in collapse for 
their most extreme combination of DM parameters. Instead, 
the excess heat goes mainly into dissociating H2, and the 
continuum-dominated cooling regime begins earlier. Never- 
theless, their work was a one-dimensional model, and fu- 
ture work should confirm that their results would hold in a 
three-dimensional study. For our current study, however, we 
do not directly follow the gas and protostellar evolution on 
these sub-sink scales, and instead aim to constrain the total 
DM reservoir available within the sink region of 50 AU over 
the initial accretion phases of the Pop III protostars. With 
these uncertainties in mind, our set-up also allows us to find 
an upper limit to the mass reached by the Pop III stars 
under various DM conditions, since we do not include the 
feedback effects of protostellar radiation on the surrounding 
gas (e.g. Smith et al. 2011; Stacy et al. 2011b). 



2.2 DM Initialization 

2.2.1 Density Initialization 

The DM, w hich was unresolve d on sink particle or stellar 
disk scales in IStacv et aLl(|2010l ). was initialized with two dif- 
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Figure 1. Density of gas and DM as a function of distance from most massive sink, as measured from the initial cosmological simulation. 
Recall that 1 GeV cm"'' corresponds to 1.7 X 10"'^* g cm~^. Blue dots represent individual gas particles, which roughly follow an 
density profile (dashed line). Solid line denotes the radially averaged DM profile to which the more refined small-scale simulations were 
normalized. 




Figure 2. Left: Enclosed mass A/cnc versus radius after refinement for Simulation A (puM r ^). Right: Enclosed mass versus radius 
for Simulation B (pdm r^^). Rod lines show the A/onc profile to which the DM was initialized. 



ferent density powerlaws , pdm oc r^^ ('Simulation A') and 
Pdm oc r^^ ('Simulation B'), with the normalization deter- 
mined by the total mass of DM found in the corresponding 
region of the original simulation. The corresponding profiles 
of enclosed mass M^nc can be seen in Figure O 

The density profiles were generated starting from an ini- 
tial uniform density field. This field was generated by placing 
particles at glass-like positions, which was achieved by allow- 
ing randomly placed particles to evolve under an artificial 
negative gravitational force until a quasi-equilibrium con- 
figuration is reached (|Whiti[l99i). Note that this method 
avoids small-scale fluctuations in the relative distances be- 
tween particles. This is an improvement upon a Monte Carlo 



sampling of the density fleld, which would be subject to such 
Poisson noise. 

The particles of uniform density po can then be trans- 
formed to a powerlaw density p oc f~" through the coor- 
dinate transformation [r,6,cj>) — >■ (f,6,<j)). The new coordi- 
nates will satisfy 

p{f)f^smO df AO d4> = por^sin6' dr A6 A(f), (1) 
from which we can derive the relation 

r oc r^/<^""', for < n< 3. (2) 

We initially align the DM density peak with the gas 
density peak, represented by the first sink particle. Figure [T] 
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shows the initial DM density profile, taken from the original 
full-scale simulation, used to normalize the DM mass in our 
current calculations. We use density units of 1 GeV cm~'^, 
which corresponds to 1.7 x 10~^* g cm"''. 

We note that our DM initialization also accounts for 
gravitational contraction on scales greater than ~ 10* AU 
(0.05 pc) during earlier phases of the central gas and DM 
evolution. This is because we normalize the DM mass based 
on the central DM content built up in the full-scale simu- 
lation only after the gas has collapsed to its maximum re- 
solvable density. Further gravitational contraction on scales 
less than lO** AU is not resolved, however, so Simulations A 
and B are initialized to cover a range of possible DM con- 
figurations on small scales in order to follow the continued 
gravitational accretion of DM onto the protostellar regions. 

For Simulation A, we use 128^ DM particles, and the 
mass of each DM particle is slightly smaller than the gas, 
rriuM = 0.01 M0. We thus normalized the DM density pro- 
file to provide a total of ~ 2 x lO** M© of DM enclosed within 
a 20 pc cube, the same DM mass found in the correspond- 
ing region of the original cosmological simulation. On these 
scales, the DM thus dominates over the 4000 Mq of gas en- 
closed within the same region. We also impose a minimum 
softening length of 10 AU. For Simulation B, we used 256^ 
DM particles, each having a mass of mDM = 1 X lO"'^ Mq. 
This yielded a total DM mass of ~ 200 Mq inside a 1 pc cu- 
bical region, the same DM mass that resided in the central 

1 pc of the original cosmological simulation. On scales of 1 
pc it is the total mass of gas, 540 Mq , that dominates over 
the DM. However, were the Simulation B profile extended 
to encompass a larger 10 pc cube, the total enclosed mass 
would be a factor of 100 greater, ~ 2 x lO" Mq, and DM 
would dominate the gas on larger scales. In comparison, the 
same volume in Simulation A holds a similar mass of DM, 
~ 1 X 10" Mq. For Simulation B we use a minimum softening 
length of 25 AU. 

A larger particle number was required for this shallower 
profile so that the DM could be resolved on scales less than 
50 AU while still extending out to distances greater than 
the stellar disk, which has a radius on the order of 1000 AU. 
Note that for Simulation B, where /9dm oc , the enclosed 
mass increases as Mcnc ~ 'TiDMA'"part oc r^, where A'^part is 
the total number of DM particles used in the simulation. We 
place approximately 15 particles within the central 100 AU 
of Simulation B (Fig. [2]). Since A'pait oc r^, this means that 
a box size that is 1000 times larger, 10^ AU, will require 
a factor of 10® more particles, A'^part ~ 1-5 x 10^. This is 
already close to the 256'^ particles available for Simulation 
B, so the half-length of the box in this case does not extend 
beyond lO'^ AU (= 0.5 pc), or a total length of 2 x 10^ AU 
{— 1 pc). For this reason Simulation B uses a box size of 1 
pc. 

For the pdm oc r"'^ profile of Simulation A, however, 
we have A^part oc r. In this case, the inner region is better- 
resolved because we now have ~ 30 particles within only 20 
AU. To resolve scales 10^ times larger than this now requires 
a factor of 10* more particles, or 3 x 10^ particles within 

2 X lO'' AU (= 10 pc). A box with this half-length will then 
use all of the available 128'^ particles, and for Simulation A 
the full length of the box is thus 4 x lO" AU (= 20 pc). 

We test multiple DM profiles since, on these very small 
scales (< 1 pc), the DM density is still a matter of un- 



certainty. Many numerical studies have predicted that an 
inner DM cusp will form within galactic haloes (e.g. 
iNavarro et al.l 1 19961 ). and that the total star-|-DM mass 
within the central regions will remai n stable to subsequent 
star formation and mergers (e.g., iLoeb fc Feeble^ l2003l : 
iGao et al]|2004l ). However, analytic and numerical models of 
DM evolution during the cooling and infall of baryons find 
that adiabatic contr action will further incre a se the central 
DM densities (e.g., iBlumenthal et al] Il986l : iGnedin et al.l 
l2004h . This is in contrast to observations of galaxies that 
point to central DM concentrations that are flat instead 
of peaked (the widely-known 'c usp-vs.-core' problem, e.g. 
pMoorg 1994; Bu rkcrt 1995; do Bl ok et al.l200ll:lGentile et al.l 
I2OO5I : ISpekkens et all [20051 : iBattaglia et all bOOSl ). The in- 
nermost resolved DM of our initial cosmological simulation 
in fact shows a proflle of approximately pdm oc r^^'^ (see 
Fig. [T]). Our simulations therefore cover a range of these 
possibilities. 



2.2.2 Velocity Initialization 

For a given spherically symmetric density proflle, an 
isotropic distribution function (DF) for the DM par- 
ticles can be generate d using Eddington's formula 
iBinnev fc Tremainjliooj . Equation 4.46): 



d*' dV 1 / dp 
- *' d*^ ^ 7f Vd*' ;*'=(, 



(3) 



where f{£) is the distribution function in units of mass per 
phase space volume, e.g., g cm""' (cm s~^)~"', and 'I' is the 
relative potential, which can be set to the negative of the 
gravitational potential as measured from the edge of the 
system, p is the mass density of the system at the given 
vp, and £ is the relative energy of DM per unit mass. If ^ 
and the relative energy £ are known for a particle, then its 
velocity v can be calculated using 



(4) 



For the input density we assumed the DM profile is cut 
off at an outer radius of 10 pc for Simulation A and 0.5 
pc for Simulation B. We further assume the profile flattens 
to a uniform-density core in the central 5 AU, well inside 
the extent of the sink particles. To assign velocities to each 
particle, we flrst divided the DM proflle into 3000 radial bins, 
centered upon the densest particle. We then calculated the 
relative potential ^ at each bin. Each particle was assigned 
the value for ^ corresponding to its bin, and a value for its 
relative energy £ was randomly drawn from the DF. The 
amplitude of the velocity v was then found using Equ. 4. 
In the Appendix we provide more details about randomly 
drawing a value of £ from the distribution function f{£) . 

To determine the velocity component along each Carte- 



sian axis {v^ 



and Vz), we next picked two angles at 



random, 6 and <j) (e.g. IWidrowl [2OO0I ). Given two random 
numbers p and q ranging between and 1, we set — 
cos"^(l — 2p) and (j) = 2-Kq We used these angles to map 
each particle in velocity space, where Vx = v sm9 cos(j), 
Vy = V sinO sincf), and Vz = v cosd. The resulting radial veloc- 
ity distribution can be seen in Figures [3] and |4] 
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Figure 3. Velocity distribution of DM at the end of Simulation A. Solid line is for particles inside 10* AU from the central DM peak. 
Dashed line is for particles greater than 10* AU from the peak. For comparison, the corresponding red lines show the velocity distribution 
to which the DM was initialized. Left: 'No-gas' case. Right: 'With-gas' case. While the velocity distribution is stable for greater than 
5x10* yr (~ 1000 free-fall times for typical peak DM density poM = 10^^ GeV cm~'^) in the 'no-gas' case, the interaction between the 
central DM and the stellar disk in the 'with-gas' case has greatly flattened the velocity distribution within 10* AU, while little change 
occurs beyond the disk at distances greater than 10* AU. 
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Figure 4. Velocity distribution of DM at the end of Simulation B. Notation is the same as in previous figure. Left: 'No-gas' case. 
Right: 'With-gas' case. The velocity distribution is stable for greater than 5x10^ yr (~ 100 free-fall times for typical peak DM density 
Pdm = 10* GeV cm"'') in 'no-gas' case. Similar to Simulation A, the interaction between the DM and the stellar disk in the 'with-gas' 
case has greatly flattened the inner velocity distribution. 



3 RESULTS 



3.1 DM Evolution 



We describe results for two sets of simulations: 'Simulation 
A', in which the DM was initialized with pum oc r~^, and 
'Simulation B', which was initialized with pdm oc r~^. 

We first describe Simulation A. We initially evolve Sim- 
ulation A with only DM for ~ 50,000 yr (~ 1000 free-fall 
times for typical peak DM density pdm = 10^^ GeV cm ) 
. The resulting DM configuration is then used as the initial 
conditions for the 'with-gas' case, in which the gas and ini- 
tial sink are placed in the center of the DM minihalo and 
evolved for a total sink accretion time of tacc ~ 20, 000 yr. 
In the 'with-gas' case, several more sinks form and a disk ~ 



1000 AU develops. For comparison, we also evolve the same 
halo for another 20,000 yr without gas (the 'no-gas' case). 

Figure|3]shows how the distribution of the magnitude of 
DM radial velocity evolves for the 'no-gas' case compared to 
the 'with-gas' case for Simulation A, while the correspond- 
ing density evolution can be seen in Fig. [S] In the 'no-gas' 
case, the velocity distribution remains stable, as shown by 
the overlap of the final velocity distribution (black lines in 
Fig. [3| with the initial velocity distribution (red lines in Fig. 
[Sjl. While the velocity distribution is very mildly dependent 
upon radius, the inner and outer velocity distributions are 
still very similar to each other since the density is nearly 
an isothermal profile. The density structure is correspond- 
ingly stable, remaining in an oc profile while forming a 
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small ~ 5 Mq core corresponding to the size of the minimum 
softening length. 

This stability of the density distribution is also shown 
in Figure [6l which displays the mass enclosed within a range 
of radii at various times in the simulation. While the 'no- 
gas' case has some initial mass increase in the innermost 
regions as the core forms, it remains stable thereafter. This 
outcome is expected, since the minihalo in the 'no-gas' case 
was initialized to be a in a stable configuration. In contrast, 
the 'with-gas' case of Simulation A shows significant evolu- 
tion due to scattering of the central DM by the star-disk 
system. The velocity distribution fiattens (Fig. [3)), and the 
central core has slightly less mass than the corresponding 
region of the 'no-gas' case, though the interaction with the 
gas is insufficient to entirely erase the small DM core. Fig.[S] 
displays how this decrease of enclosed DM mass also occurs 
on larger scales of a few thousand AU after several thousand 
yr of interaction with the central gas. 

To determine how the DM evolution in Simulation A 
varies with resolution, we also performed the same calcu- 
lation but with a smaller minimum gravitational soften- 
ing length of 5 AU (Simulation A2), and a larger softening 
length of 50 AU (Simulation A3). We find the DM evolution 
remains largely unchanged, particularly for Simulation A2 
(Fig. [7|. In the 'no-gas' cases a small central core still forms 
while the density and velocity structures are otherwise sta- 
ble for over 50,000 yr. The size of the core scales with the 
size of the softening length, and thus was largest for Sim- 
ulation A3. In each case scattering by the star-disk system 
lowers the DM density and enclosed mass. However, this ef- 
fect is somewhat different for Simulation A3, in which the 
large softening length allowed the small ~ 10 Mq core to be 
completely wiped out by the star-disk system. In contrast, 
the corresponding core in Simulations A and A2 remained 
intact but unaligned with any of the sinks. 

On larger scales between 10"^ and 10^ AU, the density 
and enclosed mass in Simulation A and A2 is actually lower 
than that in Simulation A3. Of particular interest is the 
higher ambient DM density of 10^ GeV cm"^ within 10'^ 
AU at the end of Simulation A3. This is largely because 
the accretion history of the main sink in Simulation A3 is 
distinct in that, for the latter half of the simulation, only 
a single ~ 30 Mq sink persisted while the other sinks were 
lost to mergers. The lack of multiple sinks weakened the 
rate at which the surrounding DM was scattered (see Fig. 
O, though the disk gas alone managed to significantly lower 
the surrounding DM density. Nevertheless, we find that our 
results converge for softening lengths of 10 AU or smaller, 
and even with a larger softening length the qualitative effects 
of DM on the gas and protostellar evolution will still remain 
unchanged, as will be discussed in Sections 3.3 and 3.4. 

The maximum DM density in Simulation B is much 
lower, so we initially evolve the DM-only set-up for a longer 
time of 5x 10^ yr (~ 100 free-fall times for typical peak DM 
density pdm = 10* GeV cm~'^) so that the central DM may 
still evolve for many free-fall times. We then add the gas to 
the resulting DM configuration and again evolve the simu- 
lation for face ~ 20, 000 yr ('with-gas' case), and also evolve 
the DM-only case for another 20,000 yr for comparison ('no- 
gas' case). As shown in Figure [H the DM velocity configu- 
ration remains stable for over 5x10^ yr in the 'no-gas' case, 
while the inner velocity distribution significantly flattens for 



sink 


*form [yr] 


A/flnal [Mq] 


n„it [AU] 


f-flnal [AU] 


1 





28 








2 


740 


17 


130 


260 


3 


5540 


1.0 


850 


170 


4 


9430 


0.5 


660 


5800 



Table 1. Formation times, final masses, distances from the main 
sink upon initial formation, and distances from the main sink at 
the final simulation output in Simulation A. We include the sinks 
still present at the end of the simulation (20,000 yr). 



the 'with-gas' case. We can also see that, for this DM density 
profile, the initial velocity distribution differs significantly 
depending upon radius. We find that gravitational contrac- 
tion in fact leads to an increase in DM density within the 
central 10* AU over the first ~ 10,000 yr, and the powerlaw 
of the central profile correspondingly steepens. However, as 
shown in Figure|B]and later in Section 3.3, this gravitational 
contraction is not sufficient to lead to significant amounts of 
DM within the sinks themselves. Scattering by the motion 
of the star-disk system prevents the central DM from main- 
taining these enhanced densities, and by 20,000 yr it in fact 
leads to a decrease again in the central DM densities (Fig. 

The enclosed DM mass within the central 1000 AU has 
decreased as well (Fig. |6]). 



3.2 Protostellar Mass Growth 

We also compare the growth of the sinks without DM re- 
finement ('unrefined') to that found in Simulation A and B 
(Figure [S|. The largest sink was from the 'unrefined' simu- 
lation, which grew to 34 Mq. This was still similar to the 
largest sink masses attained in Simulation A and B, 28 Mq. 
Furthermore, a total of four sinks remained at the end of 
Simulation A, as compared to three sinks in Simulation B 
and two sinks in the 'unrefined' simulation. In all three cases 
the protostellar system is dominated by a massive binary 
where each member has mass of ~ 20-30 Mq , and each case 
furthermore has a total sink mass of ~ 50 Mq . The presence 
of refined DM thus had minimal effect on the overall gas ac- 
cretion rate of the protostars, particularly for the first 5000 
yr (~ 100 free-fall times for sink density n — 10^^ cm~^). 
Instead, in all cases the sink accre tion rate at early times 
was similar to that found in, e.g., IStacv et al.1 (|2010l ) and 
IStacv et al.1 (|201ll ). 

The accretion histories are not all perfectly identical, 
however. The early mass growth of the largest sink in each 
simulation is largely driven by mergers with other sinks, 
as is evident by the steps in mass growth visible in Figure 
[8] Thus, even though each simulation began with identi- 
cal initializations of the central gas, the slight changes in 
the gravitational potential due to the different DM set-ups 
were sufficient to cause deviations in the sink merging and 
growth history. In particular, a merger at around 8000 yr 
causes the 'unrefined' case to diverge from Simulations A 
and B. The sink characteristics are summarized in Table [T] 
for Simulation A, Table[2]for Simulation B, and Table|3]for 
the 'unrefined' simulation. An example of the asymmetric, 
disk-like structure from which the sinks accrete is shown for 
Simulation A in Figure [Q] 
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Figure 5. Evolution of density structure of DM as measured from ttie densest DM particle. Red lines are analytic fits the to density 
profile to which the DM was initialized. Black lines represent 'no-gas' cases, while blue lines represent 'with-gas' cases. Dotted lines 
represent typical density threshold for DM scattering accretion to support a star, given m^vimp = 100 GeV (Freese et al. 2008b, Yoon 
et al. 2008, locco 2008, locco et al. 2008, Spolyar et al 2009). Left: Simulation A. The 'no-gas' case is shown at 50,000 yr (black dotted 
line), and after a further 20,000 yr (black dashed line). Also shown is the DM density for the 'with-gas' case after t^cc = 10,000 yr 
(blue dashed triple-dot line) and tacc = 20,000 yr (blue long-dashed line). Right: Simulation B at 5 X 10^ yr (black dotted line), and 
after a further 20,000 yr (black dashed line). 'With-gas' is shown at face = 10,000 yr (blue dashed triple-dot line) and tacc = 20,000 yr 
(blue long-dashed line). The central regions are stable in each 'no-gas' case. However, the corresponding regions in the 'with-gas' case of 
Simulation A has a lower mass and density due to scattering by the star-disk system over 20,000 yr. In contrast, the 'with-gas' case of 
Simulation B sees a temporary increase in central DM density, but this increase is insufficient to affect the evolution of the primordial 
gas, and by 20,000 yr pdm has decreased again. 




Figure 6. Enclosed DM mass versus radius as measured from the densest DM particle. Red lines show the analytic fits to which the 
DM was initialized. Left: 'No-gas' case of Simulation A at 50,000 yr (black dotted line), and after a further 20,000 yr (black dashed line). 
Note that these lines overlap. Also shown is the enclosed DM mass for the 'with-gas' case after tacc = 10,000 yr (blue dashed triple-dot 
line) and tacc = 20,000 yr (blue long-dashed line). While the enclosed mass within the inner regions slightly increases and stabilizes for 
the 'no-gas' case, the action of the stellar disk in the 'with-gas' case causes the enclosed mass to actually decrease in the inner regions 
with respect to the evolved 'no-gas case' (black), instead becoming relatively more concentrated at farther distances. Right 'No-gas' case 
of Simulation B at 5 X 10^ yr (black dotted line), and after a further 20,000 yr (black dashed line). Note that these lines overlap. Also 
shown in the enclosed DM mass for the 'with-gas' case at tacc = 10,000 yr (blue dashed triple-dot line) and 20,000 yr (blue long-dashed 
line). Some gravitational contraction of gas can be seen in the inner 10* AU at 10,000 yr, but this is insufficient to increase the DM 
concentration within any of the sinks, and the DM concentration decreases again by 20,000 yr. 
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r [AU] r [AU] 

Figure 7. Radial distribution of DM at the end of Simulation A2 and A3, which have softening lengths of 5 AU (dotted lines) and 50 
AU (dashed lines). Black lines are the state of the DM at the end of the 'no-gas' cases, after a total of 70,000 yr. Blue lines show the 
DM at the end of the 'with-gas' cases, which were each run for 50,000 yr with DM only, and another 20,000 yr after addition of the gas. 
Red lines show the analytical radial profiles to which the DM was initialized. Left panel shows the density profile, while the right panel 
shows the profile of enclosed mass. The initial DM core is largest for the 50 AU case. Regardless of softening length, however, scattering 
by the star-disk system shifts most of the central DM to lower densities. 



sink 


*form [yr] 


A/final [Mq] 


r-init [AU] 


J-flnal [AU] 


1 





28 








2 


330 


25 


66 


600 


3 


6200 


1.0 


64 


1500 



Table 2. Same as Table 1, but for sinks remaining at the end of 
Simulation B at 20,000 yr. 



sink 


tform [yr] 


Mfinal [M0] 


rinit [AU] 


rfinal [AU] 


1 





34 








2 


1160 


16 


140 


840 



Table 3. Same as Table 1, but for sinks remaining at the end of 
the 'unrefined' case at 20,000 yr. 



3.3 DM Annihilation Heating on Sink Scales 

3.3.1 Methodology of Heating Rate Measurement 

We now discuss how we determined the eflFect of DM annihi- 
lation on the central gas. We record the number of DM par- 
ticles within each sink at every timestep, and from this we 
can estimate the DM heating rate. This rate approximately 
is: Fdm oc pumi'^av) /mwiMP , where {aav) is the product 
of the WIMP annihilation cross section and relative veloc- 
ity, averaged over the WIMP momentum distribution, and 
JTiwiMP is the mass of a single dark matter WIMP. Fdm 
has further dependencies upon the energy spectrum of the 
photons produced by WIMP annihi lation, the scattering o f 
photons within the gas, etc. (see e.g. lNataraian et al.ll2009l ). 

We express the evolution of the DM annihilation heat- 
ing rate rDM(t) within the sinks in terms of a reference 
heating rate Fdmo, 



FDM(t) _ (p(t)'}sink ^ C{t) (p(f))Lk 

Fdmo (p(to)2)sink c(fo) <p(to))Lk ^ ' 

where the brackets indicate a spatial average on the scale 
of the sink radius. We choose the reference time to to be 
the time at which the average DM density within the main 
sink reaches its maximum, i.e., {/9(to))sink ~ max((p(t))sink). 
The DM clumping factors c{t) = (/5(i)^)sink/{p(i))sink and 
c{to) = (p(to)^)sink/(p(to))sink describe the variance of DM 
density fluctuations, both resolved and unresolved, on scales 
below the sink particle radius at the respective times. 

In our simulations we only track the evolution of the av- 
erage DM densities within the sinks, (p(t))sink, at high time 
resolution. On the other hand, the DM clumping factors are 
only known for the few time steps at which a full output 
of the simulation data is performed. We therefore estimate 
the ratio of DM heating rates using the ratio of the aver- 
age DM densities {p(i))sink/{p(*o))sink- This means that we 
implicitly ignore that DM clumps on sub-sink scales. 

However, even if DM clumps significantly on sub-sink 
scales, computing the ratio of DM heating rates using only 
the average DM densities is a good approximation if, as one 
can reasonably expect, the DM dumpiness on sub-sink scales 
is largest when the average DM density within the sink is 
largest. In this case, i.e., for c{t) ^ c(io), any decrease in 
the ratio of the squared average DM densities conservatively 
underestimates the decrease in the ratio of DM heating 
rates, because FDM(t)/FDMO < (p(i))sink/{p(*o))sink- Com- 
putations of the DM clumping factors on sub-sink scales at 
times where a full data output is available confirm the va- 
lidity of this assumption. 

3.3.2 Evolution of the Heating Rate 

From lNataraian et all (|2009l ) we find that Fdmo is of order 
10~* erg cm~''s~^. This is the heating rate they found at 
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Figure 8. Growth of the main sink over time in the various 
simulations. Solid line represents the simulation in which the DM 
was not refined. Dashed line is taken from Simulation A, and 
dotted line is the result for Simulation B. The final masses reached 
by the largest sink is similar in each of the cases, though a merger 
at 8000 yr causes a somewhat larger final mass for the 'unrefined' 
case. 



20 AU (10"* pc) from the central peak, given the typical 
inner density profiles of their 'fit #2' where pom — W^^ 
GeV cm~^ at this distance, and assuming a WIMP mass of 
w-wiMP = 100 GeV. These inner density profiles are simi- 
lar to the initial density profile used in our Simulation A, 
and we can thus assume a similar heating rate at 20 AU 
from the DM peak, wh ich initially lies within our main sink. 
iNataraian et al.l (|2009l ) find the DM heating rate will exceed 
the gas cooling rate within these inner regions, and due to 
the similar DM construction the same applies to sub-sink 
gas in Simulation A. 

As the DM density profile evolves, however, the heating 
rate evolves as well. As can be seen in Figur e [TOl while the 
heating rate Fdmo of lNataraian et al] (|2009l 'l may apply at 
early times for our Simulation A, Fdm within race of the 
two largest sinks shows a steady decline over their accre- 
tion time as the DM central peak is scattered away. For the 
largest sink of Simulation A we see a four order of magni- 
tude drop in heating rate, from 10~* ergcm~^s~^ to 10~^^ 
erg cm~'^s^^ by 5000 yr (~ 1000 free-fall times for typical 
peak DM density pdm = 10^^ GeV cm"''), and to zero by 
20,000 yr. This same zero heating rate was measured at the 
end of Simulations A2 and A3. 

More precisely, a measured heating rate of zero is simply 
a lower limit, as the sinks could lie within regions of DM with 
densities as high as ~ 10^" GeV cm"'^ without the discrete 
DM particles of the simulation falling within the sink radius, 
thus yielding a zero measurement when the physical heating 
is small but non-zero. Even at 10^° GeV cm , however, the 
annihilation heating rate becomes negligible. Furthermore, 
the distance between the densest DM particle and main sink 
has grown to ~ 10* AU by the end of the Simulation A. 
Figure [5] shows that at these distances the DM densities are 
much lower than 10^" GeV cm~^ at 20,000 yr, and were in 
fact lower than this even at the beginning of the simulation. 
Even without any change of the DM profile from its initial 



configuration, the motion of the sinks from the DM peak was 
sufficient to severely mitigate DM heating within the sinks. 
With the flattening of the DM profile, however, DM heating 
would not become significant unless the sinks again moved 
to within ~ 100 AU of the DM peak. Such realignment is not 
seen in the simulation, however, due to the orbital motion 
of star-disk system. 

We note that for a smaller WIMP mass, e.g. mwiMP = 
10 GeV, Fdm will be a factor of ten greater. In this case 
the DM heating rate over the first 5000 yr evolves from 
approximately 10"^ ergcm~'^s~^ to 10~^^ e rgcm~'^s~^ at 20 
AU f rom the sink's center (see discussion in lNataraian et all 
l2009f ). This roughly corresponds to Fdm falling from ten 
times the gas cooling rate to 1000 times smaller, given a gas 
density of 10^^ cm"''. The cooling rate is even greater for 
denser gas closer to the center of the sink, while the ambient 
DM density and heating rate within the sink will continue 
to decline over time. 

Along with the instantaneous relative heating rate 
shown in Figure fTOl we also compare the time-averaged heat- 
ing rate FDM.avg in Figure [TT] which for each timestep was 
calculated as the average of the relative heating rate up 
to the given time in the simulation, {p(tacc))sin]j/(p(io))sinki 
where 

(p(tacc))Lk = 7^ {p{t)}L^ dt. (6) 

Because the presence of DM within the sinks declines until 
no DM is left in these regions by the end of the simulation, 
the average heating rate over the first 20,000 yr is an order 
of magnitude less than its value at the beginning of the 
simulation. A very similar trend of DM heating rates which 
rapidly decline to a measured zero value is also seen for 
both smaller and larger resolution scales (Simulations A2 
and A3). 

We also note that a few of the sinks formed later in 
Simulation A initially developed within the persistent ~ 5 
M0 DM core. This can be attributed to the gravitational 
enhancement in this region. A typical example of the evo- 
lution of the DM heating rate within such a sink is shown 
as the dash-dot line in Figure 1101 The motion of the sink 
through the disk causes a separation between the sink and 
DM peak after less than 5000 yr, and the DM heating rate 
correspondingly declines. Thus, even when a sink is formed 
in the location most conducive to DM heating and accretion, 
large DM heating rates powered by gravitational accretion 
will be very short-lived compared to the lifetime of the star. 

Even when the softening length was varied in Simula- 
tions A2 and A3, the results were quite similar. The DM 
heating rate within the main sink dropped three to four 
orders of magnitude after ~ 10,000 yr, and no DM parti- 
cles were present within the sink radius by ~ 20,000 yr. 
However, the phenomenon of small sinks forming within the 
small DM core did not occur for the large softening length of 
Simulation A3, so this may simply be a numerical effect. It 
is furthermore interesting that in Simulation A3, which had 
only one sink in the latter 10,000 yr, the motion of the disk 
gas and the single sink alone was sufficient to significantly 
decrease the infiuence of DM. Though there may still have 
been multiplicity in Simulation A3 on unresolved sub-sink 
scales, this result shows that disk formation and subsequent 
gas rotation can minimize DM effects even without stellar 
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multiplicity. However, Figures [5] and [7] show that the reduc- 
tion in DM densities between distances of 10'^ and 10^ AU 
from the DM peak is strongest in the case of a multiple 
system. 

For the largest sink of Simulation A, we thus find that 
DM heating may be important for ~ 5000 yr, and gravi- 
tational accretion of DM will end prior to 20,000 yr. After 
this point the DM peak is well over 50 AU away from the 
protostar, beyond the edge of the sink. Note that in de- 
termining this time estimate we did not account for loss 
of DM from the simulation after it entered the sink region 
and possibly underwent annihilation. The total mass of DM 
within the whole simulation box thus never decreased, keep- 
ing the central DM potential at a maximum. We estimate 
the mass loss of DM through annihilation using the above- 
quoted heating rate, 10~* erg cm~^s ~^, which was found as - 
suming {cjav) = 3 x 10"^*^ cm^ s"^ (|Nataraian et al.ll2009l ). 
If DM within a sphere of radius 20 AU were lost at this 
rate over a period of 20,000 yr, this would correspond to 
a total DM mass loss of > lO^'^ Mq. This is a negligible 
effect, even considering that the true mass loss may be sev- 
eral times higher than this since annihilation does not heat 
the gas with perfect efflciency. The WIMP annihilation cross 
section is thus unimportant compared to gravitational scat- 
tering in determining the evolution of the DM density and 
heating rate. 

For the lower DM densities of Simulation B, typical 
heating rates within the main sink should be much lower, 
and we find the heating rates should never become signif- 
icant to the thermal evolution of the gas. The initial DM 
density profile of Simulatio n B most resembles t he shal- 
lowest profile of 'fit #1' in iNataraian et al] (|2009l ). where 
Pdm — 10* GeV cm~'^ at 20 AU from the peak. In this case 
rDM(t)/rDMO does not start at a value of one because the 
maximum heating rate Fdmo is not reached until tacc ~ 1000 
yr. The sink and DM density peak are initially aligned, and 
gravitational contraction causes the DM peak to increase 
in density by another factor of ~ 10 for the first 1000 yr. 
Given the highest attained density of ~ 10^ GeV cm~^ at 
20 AU from t he peak, this correspo nds to Fdmo 10^^'' 
erg cm^^s^^ (|Nataraian et al]|2009l ). Even when the simi- 
lar DM profile is extrap olated to nearly protostellar scales in 
INataraian et al.l (|2009l ). the heating is still not great enough 
to influence the gas collapse. Increasing Fdmo a factor of ten 
by assuming a smaller WIMP mass of mwiMP = 10 GeV still 
yields a heating rate several orders of magnitude lower than 
the cooling rate of gas within the sink. 

As the gas and DM evolve in Simulation B, this remains 
the case. The right panel of Figure [TD] shows that the heat- 
ing rate within the two largest sinks does not grow after 
tacc ~ 1000 yr. There is no further increase through, e.g., 
gravitational contraction. The growth of {p(t))sink instead 
halts, and is later reversed. By 5000 yr FDM.avg/FDMO be- 
gins a steadily decline fFig. llip . Similar to Simulation A, the 
distance between the densest DM particle and main sink has 
grown to ~ 10* AU by the end of the simulation. Thus, for 
this particular DM configuration, we find DM annihilation 
heating will not play a significant early role in the forma- 
tion of massive Pop III stars. We next discuss the possible 
longer-term effects of DM scattering accretion. 



3.4 DM Scattering Accretion on Stellar Scales 

The rate at which DM will be captured by the star through 
scattering is roughly pr oportional to the DM density in the 
vicinity of the star ('e.g. lloccd[20o3 : llocco et al.ll2008l ). Thus, 
the DM scattering rate and capture rate will evolve with 
Pdm, so we roughly calculate the evolution of the capture 
rate C as C oc (p(t))sink/(p(to))sink (blue lines in Fig. [TO)) . 
After 5000-10,000 yr the capture rate drops to 10"^ - 10"^ of 
the maximum value, and for both cases the measured value 
for C drops to zero by 20,000 yr. We again note this is a 
lower limit and that C may simply be small but non-zero. 
We give further estimates for the true physical values of C 
below. 

While the capture rate will be greatly reduced as the 
stars move to lower-density regions of DM, we also note that 
this rate depends on the velocity of the stars themselves 
with respect to the DM halo. This de pendence on veloc- 
ity lms_been_calcula^ed bv iGouldl (|l987t l and also discussed 
in llocco et all (|2008h . though we note a sepa r ate s tudy of 
scattering accretion bv ISivertsson fc Gondolol (|201ll ') which 
accounts for the fact that the concentrated DM is gravita- 
tionally bound to the star and assumes that the star does 
no t move t hroug h the DM halo. The equation as determined 
bv lGouldl (|l987i ) is: 



C^Atx [ dRR- 
Jo 

where 



dC{R) 
dV ' 



(7) 



dV V TT / M„ mwiMP 



and F is defined as 



F 



^— { (^A+A. - 1) [x{-V,v) - X{A-,A+)] 



2vA^ {\^+'^- 2 
+ '-A^e'---lA.e--l-,e--'} 



and 



3v'^{R)lJ. 



A±=A±ri, V=\l7^, 



(9) 



(10) 



X{a,b) = 



dye-" = y^[eri{b) - eTi{a)] 



R, is the stellar radius, ao is the WIMP-proton elastic scat- 
tering cross section, A„ the atomic number of stellar nuclei, 
AI„ the atomic mass of stellar nuclei, v the WIMP veloc- 
ity dispersion, v{R) the escape velocity at a given radius R 
inside the star, jj, — mwiMp/M„, and fi- — {fi ~ l)/2. As 
discussed in, e.g., Freese et al. 2008b, for Vt = 10 km s^^, 
mwiMP = 100 GeV, and v = 10 km s~^, the stellar velocity- 
dependent part of the capture rate (the factor F) amounts 
to a value of order one. For example, given a 30 Mq star that 
has contracted to 5 Rq , if the star's motion through the disk 
leads to a change in from 1 km s~^ to 10 km s~^, then F 
increases from 0.24 to 1.6 This corresponds to an increase in 
the capture rate by a factor of ~ 7 due to increased velocity 
of the star through the DM halo, llocco et al.1 (|2008l ) provide 
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Figure 9. Gas density structure of the central 10,000 AU of Simulation A after tacc = 20,000 yr. The main sink is denoted by the 
asterisk. Other sinks are denoted by diamonds. This star-disk system gradually scatters the central DM to progressively lower densities. 




Figure 10. The relative DM heating rate, approximated to be (p(t))g;j-jij/(p(to))gi„ij (black lines), within the two largest sinks over the 
20,000 yr of accretion followed in the simulation. The DM capture rate, taken as (p(t))sink/(p(io))sink (blue dotted lines), is also shown 
for the largest sink. Relative heating and capture rates are measured with respect to the highest rate attained within the largest sink. 
Left: Simulation A (Fomo ~ 10^* erg cm~^s~^). Right: Simulation B (Fdmo 10^^** erg cm^^s"^). Solid lines denote the heating 
rate for the most massive sink, the dashed lines for second-most massive sink, and dashed-dot line is for a representative small (~1 Mq) 
sink. While there is some oscillation of the DM density and heating rate within the main sink, there is a steadily decline as the star-disk 
system evolves and disrupts the initial high DM concentration. Fdm within the second largest sink stays steadily small in Simulation A. 



a simple estimate for the above capture rate given fiat cen- 
tral DM and stellar density profiles, and also assuming that 
Vf = V — 10 km s~^, values which apply well to Simulation 
A. This leads to their simplified expression: 



C 



9.2 X lO^^s"^ 



Ml 
R, 



Pdm 



10"GeV cm-» 



/ mwiMP \ ' 
lO-ascmV VlOOGeVy 



(11) 



An initial estimate for the scattering accretion rate can be 
made with the above equation, in which M« is the protostel- 
lar mass in solar masses, and R* is the protostellar radius in 



= 30 Mq 

-3 



i?* = 7 X 10^^ cm (100 Rg 



cm. We estimate A/* 

Pdm = 10^* GeV cm""^, and ttiwimp = 100 GeV. For ao 
we use the spin-dependent cross section of 10~'^*cm^, con- 
sistent with exp erimental upper lim its f or 100 GeV WIMP S 
as found by, e.g.. lDesai et al.l (|2004l ') and I Angle et al.l(l2008h, 
though we note that more recent results from lBehnke et al] 
(|201ll ') find slightly higher upper limits of ~ 7 x 10"^® cm^. 
From this we find a capture rate of 10*^ s~^. For 100 GeV 
DM particles, this corresponds to a total captured DM mass 
of 10"^ Mq over 5000 yr, and a luminosity from DM an- 
nihilation of Ldm ~ 2 x 10"° erg s"^ (= 5 x 10** L©). 
A large but brief luminosity , similar to that described in 
ISivertsson fc Gondold (j201lf l despite their different method 
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Figure 11. Evolution of the time-averaged relative heating rate within the sinks, (p(tacc))^;jj]j/(p{to))gi^]5, where (p(io))gijjij is the square 
of the maximum average density attained by the largest sink of the simulation. Left: is for Simulation A. Right: Simulation B. Solid lines 
denote the most massive sink, the dashed lines denote the second-most massive sink, and dashed-dot line is for a representative small 
{~1 Mq) sink. The sink and DM density peak arc initially aligned. At early times in Simulation B, gravitational contraction causes a 
further increase in peak DM density by a factor of ~ 10. This increase is later reversed through gravitational scattering. As the amount 
of sub-sink DM particles declines towards zero in both simulations, (^(iacc))^;^-^^. declines to an order of magnitude below (p(*o))si,j]^- 



of calculation, may then be possible. We again point out, 
however, that the initial DM capture rate drops to < 10~^ 
of the initial value by ~ 10,000 yr. This estimate also shows 
that, because DM is expected to account for only a small 
fraction of the protostellar mass itself, this effect was safely 
neglected when following the mass growth of the sinks. 

Over longer periods of time, the stars may continue ac- 
creting DM at very reduced rates throughout their main 
sequence. This includes the small, ~ 1 M© , stars that would 
need less DM to be supported by DM annihilation. There 
is a ~ 10* AU DM core of lO'^ GeV cm~^ at the end of 
Simulation A, and a similar core at the end of Simulation 
B. If this core is stable, this would yield a long-term scatter- 
ing accretion rate of lO''"^ s~^ for main-sequence stars with 
Af» = 1 M© and J?, = 1 Rq, or a luminosity from DM 
annihilation of Ldm ~ 2 x 10^^ erg s~^ (= 5 x 10~^ Lq), in- 
sufficient to power the star. Larger stars with M, = 30 M© 
and R, = 5 Rq will have a long-term scattering accretion 
rate of 2 x 10^^ s~^, and a luminosity of Ldm ~ 4 x lO^** erg 
s~^ (= 10 Lq), again insufficient to power a massive star. 
For the higher ambient DM density of lO'"* GeV cm^"^ at the 
end of Simulation A3, which had a larger 50 AU softening 
length, this would lead to a luminosity of 10-^ Lq for the 
single 30 M© star in this simulation, still a negligible frac- 
tion of the typical main sequence lumino sity of such a star 
(~ 5 X 10^ Lq, e.g. iHosokawa et al.|[2010l 'l. 

Though the precise protostellar evolution is unresolved, 
and though the contraction of the protostar may indeed 
be temporarily delayed through DM heating. Simulation A 
shows that the initial influence of DM on the protostar will 
end once the star has grown to ~ 30 Mq . Any DM that may 
have been incorporated onto the protostar early on will be 
quickly used up as an energy source for the protostar, while 
later s cattering accretion will be minimal. iNataraian et al.l 
(|2009l ) estimate the initial DM supply will be exhausted by 
~ 10^ yr, after which the protostar will subsequently un- 



dergo Kelvin-Helmholtz contraction onto the main sequence 
as usual. More precisely, they estimate the depletion of t^cp 
of DM to be 

_ Pdm ^ n-wiMP 
f-dcp — r 

Pdm PY>yi\OaV) 

lo2i\j / "IWIMP \ / PDM A""*^ 

= ^° ^^HlOOGiV J lloT^GiV^j • (^2) 

Compressing the total accreted DM mass calculated above, 
10""^ M0, to the size of the protostar (100 Rq) would yield 
a slightly higher DM density of pdm = 10^^ GeV cm . This 
would indeed correspond to a depletion timescale of ~ 10^ 
yr, significantly shorter than the expected Pop III lifetime 
of several Myr. If we determine a powerlaw fit to the latter 
10,000 yr of mass growth of Simulation A's largest sink, we 
find Af» oc t ^'^^. Extrapolating to ~ 10''' yr, we find that 
the star will have only grown to 32 Mq by the time the 
DM is depleted and the protostar can contract to the main 
sequence. Even if this depletion time were increased by a 
factor of ten such that t ~ l(f yr, the star will only grow to 
~ 40 Mq in this time, given the same extrapolation. With 
so little late-time DM scattering accretion, DM annihilation 
of the small amount of mass gathered early on is therefore 
unlikely to extend the gas accretion time of the typical Pop 
III star long enough to significantly alter the final stellar 
mass. 

While the stars should spend most of their time in the 
large low-density DM core, we have thus far neglected the 
possibility that the sinks will re-align with the smaller, high- 
density DM core (see Fig. [5]) at some points during their 
orbit through the disk. Were a ~ 100 AU region of 10^" 
GeV cm~^ to persist, a star going at an average relative 
speed of 10 km s~^ would pass through the dense region in 
50 yr, very briefly enhancing the luminosity up to 50 Lq 
and lO** Lq for a 1 Mq and 30 Mq star, respectively. This 
is less than one-tenth of the main-sequence luminosity of 
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the most massive stars in our simulations, but is much more 
significant for stars of ~ 1 Mq. However, considering the 
stellar orbits through the disk typically cover a radius of 
1000 AU, the area covered by > 10^° GeV cm"^ DM relative 
to the stellar orbital area is a factor of 10~^, so periods 
of high DM luminosity should not represent a significant 
portion of the stellar lifetime. 



4 DISCUSSION AND CONCLUSION 

We performed multiple N-body and SPH simulations with 
GADGET2 to examine the mutual interaction between Pop 
III stellar disks and central DM cores. Our results are high- 
lighted in Figure [To] We find that the motion of the Pop III 
star-disk system will cause a fiattening of the central DM 
peak. This also causes a separation between the location of 
the sinks and the DM peak of ~ 10* AU , so that after a few 
thousand years DM will not be dense enough to have an ef- 
fect on the evolution the gas and protostars within the sinks 
(Figs. [S] and [TO}. For the steeper profile of Simulation A, we 
find that even if the DM density profile were unaltered from 
its initial set-up, the relative distance between the sinks and 
the DM peak is large enough that the DM will be too diffuse 
to affect the sub-sink gas and stellar evolution. The observed 
flattening of the profile, however, allows the sinks to be as 
close as < 100 AU from the DM peak while still avoiding 
significant DM infiuence, though such close distances are 
not seen in the latter part of the simulation. DM is even 
less influential in the shallower profile of Simulation B. Our 
results show that Pop III stars will grow and contract to 
the main sequence largely unaffected by DM, and the final 
stellar masses will most likely reach only a few tens of solar 
masses. 

Previous studies of Pop III dark stars have not con- 
sidered their formation in the context of a disk and stel- 
lar multiple system, and none of these studies have utilized 
three-dimen sional calcula ti ons. O ther one-dimension a l stud - 
ies by, e.g., iFreese et al] (|2008l ) and ISpolvar et al] (|2009l ) 
found that, given a stationary star in the center of a mini- 
halo, gravitational accretion of DM will be prolonged, pro- 
viding a sufficient power source to delay protostellar collapse 
for ~ 10^ yr. This is in contrast to the abbreviated period 
of DM infiuence fo und in our s c enario of a stellar multiple 
system (< 10* yr). Ilocco et al.l (|2008l '). on the other hand, 
find that collapse will be delayed by gravitational accretion 
of DM only for ~ 10^ - 10* yr, more consistent with the 
timescales found in our work, though their calculation did 
not account for continued gas infall onto the star. 

Perhaps even more crucial, however, is the rate and pro- 
longation of scattering accretion of DM after the protostar 
has begun contraction to the main sequence. Previous stud- 
ies (e . g. Freese et al. 200 8b, lYoon et al.] l2008l : Ilocco et al] 
l2008l : ISpolvar et al.l |2009| ) have suggested that dark stars 
may survive indefinitely through scattering accretion if the 
surrounding DM medium remains sufficiently dense, usually 
PDM > 10'" - 10" GeV cm . These studies were unable to 
address how long such high densities could survive, however, 
as they could not address the evolution of the DM structure 
under halo mergers, further star formation, etc. We find, 
however, that within a Pop III multiple star-disk system 
the ambient DM density will rapidly decline, such that pdm 



becomes too low to sustain Pop III stars after only a few 
thousand yr (Fig. [5] ) . 

We note that we did not include a DM heating term 
in our simulations, and we initialized the simulations only 
after the gas had already collapsed to high density. Thus, we 
did not address how DM heating affects the initial baryon 
collapse before reaching sink densities. However, this has 
also been addressed with one-dimensional calculations (e.g. 
Spolyar et al.ll200^ : iNataraian et al.ll2009l : iRipamonti et al.l 



2010l ). These studies all found that DM annihilation heating 



will eventually dominate cooling at high densities (n ~ 10 
cm""^), scales unresolved in our study. Earlier work assumed 
that this would halt the protostellar contraction and extend 
the protostellar accr etion time. However, th e most recent 
and detailed study bv lRipamonti et al.l (|2010t ) finds that the 
collapse and thermal history of the gas will be relatively un- 
affected even under the dominance of DM heating. Neverthe- 
less, performing similar calculations in three-dimensions is 
an important task for future work. It will furthermore be im- 
portant to self-consistently follow the contraction of DM to 
small, nearly protostellar scales, particularly to better con- 
strain the inner DM density profile. We could only parame- 
terize this in the current calculation. The r esults we present 
here, however, demonstrate that even if IRipamonti et al] 
(|2010l ) have underestimated the typical effect of DM on the 
collapse of gas to high densities, the multiplicity of Pop HI 
stars leads to yet another way that the influence of DM will 
be rapidly mitigated. 

Our results are also consistent with constraints to the 
reionization history provided by WMAP and observations of 
the Gunn-Peterson trough, since a population of highly mas- 
sive dark stars (~ 1000 Mq ) only fits these c onstraints under 
a par ticular scenario of d ouble-reionization (ISchleicher et alj 
l2009t ). On the other hand. lScott et al] (120111 ) find that highly 
massive dark stars can be made consistent with our knowl- 
edge of reionization history by modifying other parameters 
of reionization models such as the star formation efficiency 
within haloes and t he escape fraction of io nizin g radiation. 
As dis cussed in, e.g.. lSchleicher et alj (|2009l ) and lYuan et alj 
(|201lh . a lack of such a dark star population also relaxes 
constraints on models and rates of DM annihilation set by 
the x-ray, 7-ray, and neutrino observatio ns in the cosmic 
background and i n the Milky Way (e . g. lUllio et alj|2002l: 



Knodlseder et al 

20071 : iMack et alibOOSl ) 



20031: iBeacom et all I2OO7I : lYuksel et al] 



In addition, our res ults easily fall in line w ith further 
constraints discussed in IZackrisson et al] l|2010h , who find 
that current survey data, including that from the Hubble 
Space Telescope, already reveals that lO'' Mq dark stars 
at 2; ~ 10 must be exceedingly rare and also unlikely to 
be observed by the James Webb Space Telescope. However, 
lower-mass dark stars are still not ruled out by observa- 
tions, and JWST may be able to detect these through the 
aid of gravitational len sing by foreground galaxy clusters 
(|Zackrisson et alj|2010l '). 

The disruption of central DM cusps on larger scales 
by baryonic mot ion and stellar feedback has been studied 
extensively (e.g. lEl-Zaiit et alj|200ll: iGnedin fc Zhao"200: 
Mashchenko et al. 20061 : Governato et al] 2010l : iCole ct al 



I 



20111) . Earlier work has found that the transfer of energy 



from the baryons to DM through gravitational interaction 
can indeed fiatten the central DM density profile, and our 
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work shows similar results on stellar disk scales. In short, we 
expect that previous studies of Pop III star formation lost 
little accuracy due to exclusion of DM annihilation effects 
at high densities. We furthermore do not expect to detect 
either low-redshift Pop III stars whose lifetimes have been 
extended through accretion of DM or Pop III stars that have 
grown to extremely large masses due to a gas accretion phase 
that was extended through DM annihilation effects. This is 
because the influence of DM on Pop III stars will already 
end well before the stars have reached more than a few tens 
of solar masses. 
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APPENDIX A: DRAWING FROM THE 
DISTRIBUTION FUNCTION 

We here describe in more detail how we randomly draw a 
value of relative energy £ from the distribution function 
f{£)- For a given density profile p(r) we generate a func- 
tion '^{r), where r is the distance from the center of the 
spherically symmetric DM structure. We set 



*(r) 



-$(r) -I- $0, 



(13) 



where $ is the gravitational potential and $o is a normaliz- 
ing constant. We determine $(r) over the range < r < _R, 
where R is the radius of the entire DM structure, and 



$(r-) = -47rG ^ / p(r')r'^dr'. 



We also set 



$0 = —AnG 



R 



p{r')r'^ dr' 



(14) 



(15) 



so that ^'(r = R) is zero. Once is known for the rele- 
vant range of radii and densities, this can be plugged into 
the expression for ,f{£). We compute f{£) over 3000 loga- 
rithmically spaced values of £, from Q ^ £ ^ 'I' max, where 

*max = -$0. 

When assigning velocities to each particle, we first de- 
termine the value for based upon the particle's radial 
distance r from the center of the DM structure. For the 
range of £ values we furthermore determine the probability 
P{£ ^ £o) that the DM particle will have a relative energy 
£ ^ £o for the gi ven ^. 

Recall from iBinnev fc Tremaind (|2008h that we are 
defining the distribution function such that 



(16) 



Because we have a sp herically symmetric system, and given 
that V = i/2(* -£), we can also write 



(17) 



= 4n / f{£')^2{^-£')d£'. 
Jo 

From this we derive the following expression for deter- 
mining P{£ ^ £o): 

P{£ < £o) = ^ f{£')^2{<f-£')d£'. (18) 

P(£ ^ £o) is only considered over the range ^ £ ^ ^I", 
such that P(£ > *) = 0. 

We furthermore normalize P{£ ^ £o) according to 

Pnorm(£ ^ £o) = P(£ =55 £:o)/iV(*), whcrC 



47r 



P(*) Jo 



/(£:') v/2(* -£')d£'. 



(19) 



The N{'i') normalization factor ensures that the total prob- 
ability over ^ £ ^ "ii will sum to one. 

With Pnorm determined, we next generate a random 
number ps between and 1. Finally, we find the value of 
relative energy £o such that Pnorm (f ^ £o) = Pe- The par- 
ticle is set to this £o, and velocity of the particle is then 
determined as described in Section 2.2. 
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